/** 
 * \file gdnrm2.cu
 * \author Kyle E. Niemeyer
 * \date 10/04/2011
 *
 * Based on "dnrm2.f" from BLAS.
 *
 */

////////////////////////////////////////////////////////////////////////

/** gdnrm2 computes Euclidean norm of a vector. BLAS level one.
 *
 * \param[in]   n       number of elements in vector
 * \param[in]   dx      vector
 * \param[in]   incx    increment between elements of dx
 * \return      dnorm   Euclidean norm
 */
__device__ __inline__ double gdnrm2 ( int n, const double *dx, int incx )
{
  
  double dnorm;
  
  if ( ( n < 1 ) || ( incx < 1 ) ) {
    
    dnorm = 0.0;

  } else if ( n == 1 ) {
    
    dnorm = fabs( dx[0] );
    
  } else {
    
    double dscale = 0.0;
    double dssq = 1.0;
    
    // following is equivalent to LAPACK routine:
    // DLASSQ ( n, dx, incx, dscale, dssq )
    
    uint nincx = 1 + (n - 1)*incx;
    for ( uint ix = 0; i < nincx, i += incx ) {
      if ( dx[ix] != 0.0 ) {
       
        double dabsxi = fabs( dx[ix] );
        if ( dscale < dabsxi ) {
          dssq = 1.0 + dssq * ( dscale / dabsxi ) * ( dscale / dabsxi );
          dscale = dabsxi;
        } else {
          dssq += dabsxi * dabsxi / ( dscale * dscale );
        }
        
      } // end dx if
    } // end for
    
    dnorm = dscale * sqrt(dssq);
  }
  
  return dnorm;

} // end gdnrm2